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A METHOD FOR REDUCING THE SENSITIVITY 
OF OPTIMAL NONLINEAR SYSTEMS 
TO PARAMETER UNCERTAINTY 

By Jarrell R. Elliott 
Langley Research Center 
and 

William F. Teague 
University of Kansas 

SUMMARY 

The parameters of a nonlinear dynamical system which is to be controlled optimally 
are not always accurately known. It may therefore be desirable to accept a reduction in 
the predicted nominal performance of the system in exchange for the ability to better pre- 
dict the outcome of the system, or plant, operation. 

In this paper relationships to predict mathematically the sensitivity of the system to 
parameter errors are derived and used to establish a procedure for reshaping the optimal 
solution to reduce the statistical uncertainty in the terminal conditions of the system due 
to known statistical characteristics of the system parameters. The procedure requires 
the introduction of an augmented performance index which is a linear combination of the 
original performance index and a positive scalar measure of the system sensitivity com- 
posed of the weighted sum of the variances of the performance index and terminal con- 
straints. The augmented performance index, the sensitivity partial derivatives, and the 
original state variables are used in the formulation of a new, higher dimensioned optimi- 
zation problem of the same form as the original problem. The new problem introduces 
certain weighting factors which permit different relative importance to be attached to 
different types of sensitivity, such as position relative to velocity, and which allow for 
adjustment of performance degradation and of sensitivity reduction. 

The procedure developed was illustrated by solving a nonlinear multiparameter 
rocket- trajectory problem. An algorithm based on the method of steepest descent was 
used to solve the problem because of the widespread use and proven versatility of this 
numerical technique. The example solutions serve to show the tradeoffs made possible 
by changing the weighting factors and to illustrate the radically different solutions one 
can obtain when sensitivity considerations are included in the problem formulation. 



INTRODUCTION 


In calculating the open- loop control time history of a physical process or plant, it is 
common to assume that the plant behavior or outcome can be predicted with the aid of a 
mathematical model of the plant with known or deterministic values of the plant param- 
eters. More often than not, however, the plant parameters are stochastic, not determin- 
istic, and serious discrepancies between predicted outcome and actual outcome may occur 
as a result of plant-parameter variations of one of two types: (1) parameter variations 
during plant operation, or (2) errors in the estimates of fixed plant parameters. These 
effects (that is, the sensitivity of a dynamical process to parameter variations) should be 
a consideration in the design of any control system. 

While many studies related to plant sensitivity have been reported in the literature 
(see refs. 1 to 14 for a non exhaustive list), most of them deal with linear problems. It 
also appears that little effort has been devoted to computational aspects (a notable excep- 
tion is ref. 14). The present paper deals with the stochastic nature of the parameters in 
nonlinear problems and shows how a well-known computational algorithm in optimization 
theory may be applied to obtain numerical results. 

The problem considered is a multivariable, multiparameter optimal control problem 
(with terminal constraints) whose mathematical model consists of a set of ordinary first- 
order nonlinear differential state equations. The analysis is limited to parameter varia- 
tions of type (2). The sensitivity of the performance index and the sensitivity of terminal 
constraints to parameter variations, that is, the partial derivatives of these quantities 
with respect to the parameters, are used to construct a scalar measure of plant sensi- 
tivity. This sensitivity measure, which is the expected value of the weighted sum of the 
squares of the sensitivity partials, is multiplied by a weighting factor and added to the 
original performance index to form an augmented performance index. Minimization of 
this augmented performance index, with appropriate weighting factors, results in solu- 
tions which are less sensitive to parameter variations. 

A simple nonlinear example problem, representative of rocket flight in a uniform 
gravitational field, is worked out in detail to show the steps required in setting up and 
solving plant sensitivity problems and to illustrate how reducing the sensitivity of a plant 
to parameter variations can modify the open- loop control time history and state- variable 
time history of the plant. The example also serves the purpose of demonstrating the use 
of the cumputational algorithm in solving a typical problem. 

This research was conducted at NASA Langley Research Center, with William F. 
Teague in residence under a grant arrangement with the University of Kansas. 
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SYMBOLS 


A n x q sensitivity coefficient matrix, 9 fy 9 a 

Aj partition of matrix A (i=l,2,...,n) 

a system parameter vector with components aj,...,aq 

a mean (nominal) value of a 

Cj) aerodynamic drag coefficient 

D aerodynamic drag 

dP "length" of control step 

o 

(dP) control step -size measure 

expectation of random variable { [ ~) 

F n x n state coefficient matrix 

Fj partition of matrix F (i=l,2,...,n) 

f l,f 2?f 3 components of _f 

f_ governing equations of state 

^ governing equations of augmented state 

g gravitational acceleration 

gj vector elements of augmented -state differential equations (i=l,2,...,n) 

1 0 ^,,1^0,1^ constant matrices used in steepest -descent algorithm and defined by 

equations (A 16) 

augmented performance index 
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modified drag coefficient 

vector set of terminal constraint values (in the example problem, Lj and 
L2 are the values for altitude and vertical velocity component, respectively) 

number of terminal constraints 

dimension of control vector u; vehicle mass 

initial vehicle mass 

time rate of change of mass 

dimension of state vector x 

number of augmented -state variables 

covariance matrix 

dimension of system parameter vector 

sensitivity matrix, 9x jd a; frontal surface area 

partition of matrix S (i=l,2,...,n) 

vehicle thrust 

time (independent variable) 

final time 

initial time 

horizontal velocity component, X3 

m -dimensional control vector 

total velocity magnitude 

vertical velocity component, X4 

m x m (symmetric) control weighting matrix 

relative weighting factors (i=l,2,...,Z+l) 



w g sensitivity weighting factor 

x n-dimensional state vector 

x ii -dimensional augmented -state vector 

X 2 >X 3 components of x 

x range (with mnnerical subscript denoting a specific state variable) 

y altitude 

y flight -path angle 

6 thrust -attitude control angle 

n-dimensional adjoint vector associated with payoff 
A (l + 1) x n variable matrix integrating factor 

l x n matrix of adjoint variables associated with terminal constraints 
i u constant Lagrange multiplier 

v_ £ -dimensional constant Lagrange multiplier vector 

p atmospheric density 

4> sensitivity performance index 

4 > p performance index (not including sensitivity) 

<p s sensitivity measure 

Z-dimensional terminal constraint vector 
i/'j components of >p (i=l,2 

( transpose of matrix ( ) 

A dot over a symbol indicates a derivative with respect to time. The symbol 6 
denotes a variation in a quantity, and A denotes a first-order perturbation. 
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PROBLEM STATEMENT 


Consider the following fixed- time open-loop optimal control problem: 

Minimize 

0p = ^p(— Of)) 
for the system 

x = f(x,H,a,t) ; x(t 0 ) =x 0 (given) 
with terminal constraints 
- L = 0 

where 

c pp is a scalar performance index, 
x is an n-dimensional state vector, 

£ is an n-dimensional state derivative vector, 
u is an m-dimensional control vector, 

a is a q-dimensional parameter vector, whose mean and covariance matrix 
are known ^E<QQ> = a ; E<(5a oa^} = p), and 

ip is an 2-dimensional terminal constraint vector. 

^All vectors are column vectors unless superscripted by T which indicates vector or 
matrix transpose, except f Rector ) ’ * s a row vec tor.) 

It is desired to combine some measure of sensitivity (p s to <pp such that when the 
combination is minimized, the performance index $p and the terminal constraints \p 
will be less sensitive to perturbations in the parameters a. It should be realized that 
this reduction in sensitivity will degrade the performance index <pp. However, since 
a system may have more performance capability than is needed, a system user may be 
willing to sacrifice some performance in order to reduce terminal perturbations or 
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errors in (pp and jp. In order to provide different degrees of performance loss and 
sensitivity reduction, a weighting constant is put on the sensitivity measure cp s , and the 
sensitivity performance index to be minimized is then defined as 

cp= <t>p + w s 4> s 

The sensitivity measure should reflect the fact that the system is stochastic. The 
sensitivity measure will therefore be taken to be an expected value of the weighted sum 
of the squares of the perturbations in ©„ and each ^ due to parameter variations. 

r 

(It is necessary to square the perturbation in cpp and in each (i=l,2,...,Z) to insure 
that the measure is positive). Now since the parameter variations are of type (2) 
described in the introduction, the perturbations to first order are 

. 9 *p „ 

= -5T 


and 


8*4 

Arp, = — - 6a 
i 8a - 

The squares of these perturbations may be expressed as 

,T 


(i=l,2, 


K) 2 = f 


and 


M 2 = 


8 *1 ..T^i 


6a 6a 


■ — oa 6a l -r— 
8a - — \ 8a 


(i=l,2,-..,0 


The sensitivity measure contains certain weighting factors, to be called relative weighting 
factors, which are required because one A^ quantity may be considered more, or less, 
important than another (for example, one may be more concerned with velocity errors 
than position errors in some application although both are used as terminal constraints). 
Finally, the sensitivity measure is 


f» s = eTw^A^p ) 2 +w 2 (A* i y 


+ w 


i+1 


(^i) 1 


E < w i-9T°» 5 a [sr 


l T 


+ w 


8*1 i 

2 ~dT 6 - 5 - 


f^l 

\ 9a 


dxp 


+ . 


I 


+ W Z+1 “Sa 


6a 6a 



( 1 ) 
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On a trajectory defined by some open-loop control time history, the partial derivatives 
inside the expected value are constant. By well-known identities involving the expected 
value (see ref. 15), equation (1) may be rewritten as 




8^1 
7 2 a a 


hi 




Since E {oa 6a "y = P, this equation becomes 


= W 


a 4> D /9 4 > 


1 aa \ aa 


SI#/, /8^A T S'f'j / 8 i M T 

* w 2 -sr p br) + --- + w !+i-§5- p br) 


where P is the covariance matrix of the random parameters which is assumed to be 
known along with the mean value, that is 


E(a} = a 

This expression for <p s (eq. (2)) may be recognized as the weighted sum of the vari- 
ances of the perturbations in performance and constraints. Finally the sensitivity per- 
formance index for minimization is 

* = - w s*s(^f) < 3 > 

Now return to the problem of computing and A»//.. Since these perturbations are 

explicitly functions of the terminal state, they may be written as 


b - llE g a - _E £= 5a 
b P 8a - 9x 0a - 


9^ fy* 0x 
Ai/q = 5a = -—=■ ■— 6a 
1 0a - 0x 9a - 


(1=1,2,...,*) 


The remaining problem is that of determining reference 16 it is shown that 

for parameter variations of type (2) and for the system differential equations 


x = = f^x,a,t), -g=(t) is the solution of 
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( 5 ) 



£l 

CD 

1 

8x 

8f 

£X/, 


IS 1 

1 

8a 

+ 8a ’ 

3a\ 


By defining 


S 


8x 

8a 


F 


8f 

8x 


and 



equation (5) may be written as 

^ S = FS + A ; S(t 0 ) = 0 (6) 

Then <£ g becomes 



A simple, easy-to-work-with form for <j> may be determined by reordering and pos- 

s 

sibly introducing some new variables such that 
4> p (x(tf)) =x 1 


and 


^i(x(tf)) =x i+1 {1=1,2,..., 1) 

This procedure may, if new variables are required, introduce a new system of differential 
equations. The notation used will remain unchanged, however, and the system of differ- 
ential equations will continue to be called x = f . 

T 

Now partition S, where (i=l,2,...,n) denotes the rows of S, as follows: 
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Then 


S = 


V 


<f> =(*l + w 1 S 1 T PS 1 + w 2 S 2 T PS 2 + ■ • • + w Z+l S Z+l TpS Z+l) 


that is, 


F = 


F^S + 

Ai T ; F 

and i 

r 1 

l 


A T " 
A 1 

F T 
F 2 


7t 

a 2 

• 

; a = 

* 

f~t 

r n 


A T 
A n 


Observe that this optimization problem involving the sensitivity index has the same form 
as the original problem involving only the performance index (but with more dimensions) 
and may be solved by any of several algorithms. Restated, the problem is to minimize 


<p = 4>(x^, s 1 ,s 2 ,. . -> s £ + ^ 

subject to the differential constraints 

x = i = i(x,u>a,t) ; x(t 0 ) =x 0 

Si = S T F 1 + - g^x^u,!^) ; Sj (t Q ) = 0 


( 8 ) 


S n = S T F n + A n = g n (x,S i ,u,a,t) ; S n (t 0 ) - 0 


( 9 ) 
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and the terminal conditions 


V^(x(tf)) = k 


( 10 ) 


The augmented-state terms are defined as follows: 



ax ^ / _ \ 

Then the differential constraints are — = f x,u,a,t) . Thus it may be seen that the only 

dt 1 

difference between this problem and the original problem is that it is now possible to 
control, at the expense of performance index <f> , the sensitivity of the trajectory, and 
the problem is now an n(l + q) -dimensional state -variable problem rather than an 
n -dimensional one. Also several weighting factors (w s and have been 

introduced into the problem. These weighting factors are to be used in applications to 
control the loss in the performance index <p^ and to properly emphasize the importance 
of one type of terminal error relative to another type. The contribution and effect of 
these weighting factors will be illustrated in an example problem. 

It should be pointed out that solving the differential equations for the sensitivity 
partials gives all the information needed to perform a first-order error analysis on any 
of the state variables of the original system. In the example problem, a comparison will 
be made between first-order error analysis results obtained by using the sensitivity 
partials and first-order error analysis results for which one-sigma errors were intro- 
duced one at a time into the equations of motion. 

The discussion herein deals with the fixed- time problem. However, fixed time was 
used for convenience and clarity only and is not a limitation of the formulation. 


While almost any of the algorithms for solving optimization problems would be 
applicable here, it was decided to use the steepest -descent algorithm as outlined in ref- 
erence 16 because of its widespread use and proven versatility. The procedure used in 
this algorithm is to linearize about a solution x (t) provided by some reasonably chosen 
control time history u*(t) (which in general neither minimizes <p nor satisfies 
- L = 0) and to solve for 6u(t) (which improves 4>(k) anc * A new solution 

provided by u*(t) + 6u(t) is obtained and the procedure is repeated successively until 
~ k is sufficiently close to zero and can no longer be decreased. This final 

solution is said to be optimal although no necessary conditions for optimality have been 
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satisfied. A brief derivation of the necessary relationships with notation common to 
many steepest-descent programs is provided in appendix A. The algorithm was applied 
to the following example problem. 

A NUMERICAL EXAMPLE 


Problem Statement 

The example problem is a fixed-time problem in which it is required to determine 
the thrust-attitude program of a single-stage rocket vehicle starting from rest and going 
to specified terminal conditions of altitude and vertical velocity which will maximize the 
final horizontal velocity. The idealizing assumptions made are the following: 

(1) A point-mass vehicle 

(2) A flat, nonrotating earth 

(3) A constant- gravity field, g = 9.8 m/sec^ (32.2 ft/sec2) 

(4) Constant thrust and mass- loss rate 

(5) A nonlifting body in a nonvarying atmosphere with a constant drag parameter 

K D =|pC D S, where S is the frontal surface area. 

The coordinate system and pertinent geometric relations and terms are shown in 
figure 1. The differential equations of motion needed in the algorithm setup are 


= — (T cos B - K n uV) = xt = f i 
dt m \ ' J t A 


dt m 

dy • f 

-_ = v =x 2 = f 2 


( 12 ) 


at = m ( T Sin 6 ” K ° UV ) ' g = *3 = f 3 


where m is the vehicle mass and where 


V 


= \/u 2 + v^ 


and 


m(t) = m 0 + mt 

An equation for dx/dt is not included in the equations of motion becuase x does 
not enter into the problem. However the equation x = u was integrated separately to 
obtain range. The equation for dm/dt is analytically integrable, since dm/dt is a 
constant, and is therefore not included in the set of differential equations. 
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The parameters of the problem, which will be considered fixed constants whose 
precise values are unknown, are thrust level, mass-loss rate, initial vehicle mass, and 
modified drag coefficient. The parameter vector is defined as 

a T = ^T,m,rn 0 ,Kq] 

These parameters are assumed to be statistically independent and to have a normal dis- 
tribution function with mean and one-sigma values as shown in table I. 

TABLE I.- SYSTEM PARAMETER VALUES 


Parameter 

Mean value 

One -sigma value, 
fraction of mean value 

Thrust, T, kN (lbf) 

17.8 (4000) 

0.0067 

Mass -loss rate, m, kg/sec (slugs/sec) 

-9.1 (-0.62) 

0.0167 

Initial mass, m Q , kg (slugs) 

1433.6 (98.259) 

0.010 

Modified drag coefficient, Kq, kg/m (slugs/ft) 

0.048 (0.001), 0.479 (0.01) 

0.0167 


These values are considered representative of the current state of knowledge in 
solid rocketry. ^The mean values used were chosen to conform to a flight of 100 seconds 
for a rocket with specific impulse of 200 seconds, an average acceleration 


1 flOO T 

100 J m(t) ^ a ra ^° to final mass m^t 0 ^m(tf) equal to e, the 

base of the natural logarithm.) 


The boundary conditions at t = 0 and t = 100 seconds are shown in table n. 


TABLE II.- BOUNDARY CONDITIONS 


Variable 

Initial conditions 
(t = 0) 

Terminal conditions 
(t = 100 seconds) 

u, m/sec (ft/sec) 

0 

Maximum 

y, m (ft) 

0 

15 240 (50 000) 

v, m/sec (ft/sec) 

0 

0 


In the notation previously introduced 
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%(k) = - u fe) 

ip 1 = y(t f ) - L x = 0 ; Lj = 15 240 m (50 000 ft) 

^2 = v (tf) - T ->2 = 0 ; ~ 0 m/sec (0 ft/sec) 

x T = [-u,y,v] 

and the vector control variable u in the general formulation is the thrust -attitude 
angle 9. This completes the problem statement (without trajectory -sensitivity 
considerations). 


Sensitivity Relations 

The sensitivity -matrix differential equation is 


— = FS + A 
dt 


where 



3u 

8u 

0u 

0U 


S T- 


" 8T 

8m 

8m o 

8K d 


°1 

s = 

9y 

_9y_ 

9 y 

9 y 


S T 
b 2 


3T 

8m 

O 

a 

CO 

8K d 



dv_ 

8v 

8v 

8v 


s"" T " 


8T 

8rti 

8m 0 

8K D J 


& 3 



r u 2 

v+ V 

0 

uv 

V 


V 

w °|s 

1! 

fe 

0 

0 

m 

k d 




uv 
' V 

0 

V + ^ 

v_ 


*7 


and 



cos 6 

tfi 

f l 

-uV 


rA T i 

a i 

A =-7k 

0 

0 

0 

0 

= 

A T 
A 2 


-sin 9 

t( f 3 + s) 

( f 3 + $ 

vV 


a" t " 

1T3 J 


(13) 


(14) 


( 15 ) 
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By using the partition notation indicated, differential equations for the row vectors may be 
written as 


dS ; 


= F* S + Aj 

dt 1 1 


or in column vector form 


dS* rp 

dT = s F i + A i 


(i=l,2,3) 


(i=l,2,3) (16) 


The sensitivity measure is 



which may be written as 
3 

4> s = 2 w i s i T PS i (t f ) 
i=l 

where 



is a diagonal covariance matrix due to the assumption of statistical independence in the 
parameters. 


Problem Statement Including Sensitivity 
The augmented -state vector is 


.t JtI.tI.tLtI 
- "If ! S 1 ! S 2 j S 3 J 




u >y,v, 


ay av 


av 


3u a'u 8u au ay 
aT’am’am 0 ’aK D ’ 9 T’- • -’aK D ’ 9 T’- • '’8 Kd| 
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with the initial conditions 


x T (t 0 ) = [6,0,. . .,§ 

The problem restatement with sensitivity considerations is to minimize 


<t> = <t>p + w s 0 s 


subject to 


dx T ! 

1 

dS^ 

idS 2 T 

|dS 3 T 

dt | 

dt 

! dt 

dt 


and 




= x 2 (tj) - Lj = 0 ; L x = 15 240 m (50 000 ft) 
ip 2 = X3^tf) - h2 = 0 ; L2 = 0 m/sec (0 ft/sec) 


ifa>) " 


= 0 


(19) 


( 20 ) 


Numerical Results 

As previously mentioned, the steepest -descent algorithm of reference 17 was imple- 
mented to obtain numerical results. These results were computed with the use of a 
fourth -order integration subroutine with a fixed step size of 0.5 second. For this example 
problem each iteration in the solution required integration of 15 augmented state differen- 
tial equations plus the forward integration of the- range equation and backwards integration 
of 45 adjoint variable differential equations. Each iteration required about 22 seconds on 
the Control Data 6600 computer system used with no special effort having been made to 
keep computer run time down. Several check solutions were computed with the use of a 
0.125-second fixed step size. These check solutions always agreed to at least five, and 
usually to at least seven, significant figures with those computed with the use of the 
0.5-second step size. 

Numerical results for a variety of cases were obtained by changing the values of 
Kj) (the modified drag coefficient), w s (the sensitivity weighting factor), and wi,W 2 ,Wg 
(the set of relative weighting factors). The different combinations for which results were 
obtained are shown in table HI. 
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TABLE ID. - CATALOG OF CASES 


Relative -weighting-factor sets 

k d 

W S 

Set (1): = wg = 1; W 2 = 10 

0 

0, 0.1, 0.396, 1.0, 10.0 


f 0 

0, 0.1, 0.3, 1.0, 10.0 

Set (2): wi = wg = 1; W 2 = 10 -4 

X 0.001 

0, 0.1, 1.0, 10.0 


The relative -weighting-factor sets were chosen on the basis of the relative impor- 
tance attached to velocity (w^ applies to horizontal velocity and W 2 applies to vertical 
velocity) and position (wg applies to altitude errors) . Set (1) gives equal weight to a 
0.305-m/sec (1-ft/sec) velocity error and a 3.05-meter (10-foot) altitude error, while 
set (2) gives equal weight to a 0.305-m/sec (1-ft/sec) velocity error and a 30.5-meter 
(100-foot) altitude error. 

The modified drag coefficient Kj^ was set at zero, simulating vacuum flight, and 
at 0.001, a representative value for small rockets in the earth's atmosphere. Setting 
Kj) at zero reduced the number of parameters to three. It also made possible, for 
w s = 0, an analytical calculus -of -variations solution for the optimal thrust -attitude time 
history; namely, the well known linear tangency law (ref. 18). 

tan 0(t) = a + /3t 


where a and /3 are constants determined by the boundary conditions of the problem. 
This calculus -of -variations solution was used to validate the steepest-descent-algorithm 
programing and solution. The comparison of results showed negligible differences. 

Also, for Kj) = 0, it can be shown that 


8y_ = I. 
8T T 



and 


S-H*)** 

Therefore, these partials will not change with w s for = 0. 

As an independent check on the sensitivity partials, perturbations in u, v, and y 
at tf due to plus and minus one-sigma errors in each parameter, taken one at a time, 
were computed. For example, the one -sigma (la) perturbation in u due to parameter 
errors was computed by using the relation 
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Ill 111!! IIIIIIHI 


Au 


( 2 ) = 1 

la 2 



M 





u 

_ 







+ 


“(tf) 


+la ' 
m 0 




+ 




1/2 


( 21 ) 


where u ( t f) + j a means 


u at on a trajectory with 


T increased by its one-sigma 


value, where u^t^ means u at t^ on a trajectory with m decreased by its one- 

~ l0 m 

sigma value, and so forth. Similar computations were made for Av^ and Ay^. One- 
sigma perturbation values computed by this method were compared with one-sigma values 
computed by using the sensitivity partial derivatives. For example, the one-sigma (la) 
perturbation in u was computed as follows by using the sensitivity partials: 


Au 


(D = 

la 


lf 5T 16 


3u * 

3rii 6m la 


8u 
8m r 


5m, 


o,la 


~ 6K n 
8K D d 


ni/2 


( 22 ) 


where ST-,-, 6m-| and so forth are one -sigma perturbations in the parameters. Simi- 
lar computations were made for Av^ and Ay^; . Superscript (1) refers to la pertur- 
bations obtained by using sensitivity partials, and superscript (2) indicates that the method 
of computation was to obtain an average perturbation value by assuming first plus and then 
minus la variations in each parameter. 


Discussion of Numerical Results 

Numerical results for the converged trajectories are summarized in figures 2 to 11, 
which give control time histories and trajectory plots, and in tables that give sensitivity 
partials, one -sigma perturbation values, and $ g and 4>p- 

Figure 2 shows the control time history used initially and the converged control 
time history after 10 iterations for the conditions indicated (Kp = 0; w s = o) to illustrate 
how a reasonably chosen control time history is modified by the algorithm to obtain an 
"optimal" solution. Optimal here means that a gradient function has been reduced several 
orders of magnitude or to some reasonable value. A plot of tan Q( t) as a function of time 
was made for this converged trajectory, and the intercept and slope of a straight line 
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fairing to that plot are compared in table IV with the values of the constants a and j3 
required in the calculus -of -variations (C.O.V.) solution. 

TABLE IV.- LINEAR -TANGENCY-LAW CONSTANTS 


Constants 

Faired values 

C.O.V. values 

Of 

3.49 

3.49 

0 

-0.0401 

-0.0404 


The agreement is considered to be good, in view of the crude fairing method and 
other factors relating to how the C.O.V. values of a and /3 were obtained. Also the 
trajectories showed good agreement. 

Table V(a) and figures 3, 4, and 5 summarize converged numerical results for 
k d = o and set (1) relative weighting factors for several values of the sensitivity 
weighting factors. The sensitivity partials of table V(a) show consistent trends, with 
altitude and horizontal-velocity sensitivity partials decreasing in magnitude while 
vertical -velocity sensitivity partials are increasing in magnitude with increases in w s . 
These trends are also reflected in the root -sum -square perturbations; one -sigma (la) 
perturbation values of altitude and horizontal velocity decrease while those of vertical 
velocity increase. The percentage of the original value (value with w s = 0) of Aui a , 
Ay-^, Avj ct , 4 > g , and is plotted in figure 3 against sensitivity weighting factor. 

For values of w s larger than about 0.3, little change takes place in Avi a , Ay 1(J , and 
<p s , but the performance continues to degrade along with a decrease in Au 1(r The 

control time histories for these cases, shown in figure 4, exhibit an interesting charac- 
teristic. As w g becomes larger, the control tends toward a bang -bang type of control 
where the thrust is either directed straight up (vertical) or straight down. While it 
appears that the steepest -descent program indicates the existence of a bang -bang optimal 
control law, attempts to predict this behavior analytically have been unsuccessful. Alti- 
tude is plotted against range in figure 5 where it may be observed that the trajectory 
becomes steeper as w s increases. 

Table V(b) and figures 6 to 9 summarize the results for Kq = 0 and set (2) rela- 
tive weighting factors. Set (2) puts less emphasis on the altitude sensitivity partials than 
set (1). Thus altitude sensitivity partials increase for set (2) rather than decrease as 
they did for set (1), and both horizontal and vertical -velocity partials decrease. Figure 6 
clearly illustrates these results in the plots of Ay^, Av.^, and Auj a - comparing 
this figure with figure 3, it may be seen that Ay 1(J . and Av Ja have switched positions 
on the plots. Also it appears that <fi s begins to level off and remain essentially con- 
stant at w c ~ 10.0 in figure 6 whereas it leveled off and became essentially constant at 
s 
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TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION 


(a) Kq =0; = Wg = 1, Wg = 10-2 (Set (1)) 


Sensitivity partials 

Values for sensitivity weighting factor w s , of — 

and related terms 
(*) 

0 

0.1 

0.396 

1.0 

10.0 

8u/9T 

1.106 

0.973 

0.372 

0.163 

0.020 

9u/9m 

-6 322 

-5 970 

-2 044 

-838 

-101 

9u/ 9m Q 

-84.98 

-77.34 

-28.07 

-11.93 

-1.44 

< 

110.1 

101.3 

36.1 

15.2 

1.8 


110.2 

101.4 

36.2 

15.3 

1.8 

9y/8T 

52.75 

52.75 

52.75 

52.75 

52.75 

9y/9m 

-96 140 

-82 950 

-73 640 

-72 790 

-72610 

9y/9m 0 

-2 755 

-2 672 

-2613 

-2 608 

-2 606 

< 

3 209 

3 100 

3 025 

3 019 

3 017 


3 210 

3 100 

3 025 

3 019 

3 017 

9v/9T 

0.805 

0.805 

0.805 

0.805 

0.805 

9v/9m 

-1419 

-1646 

-1805 

-1820 

-1823 

9v/9m 0 

-41.74 

-43.17 

-44.18 

-44.27 

-44.29 


48.57 

50.50 

51.91 

52.04 

52.07 

< 

48.57 

50.51 

51.93 

52.12 

52.18 

d> 

4 423 

3 892 

1489 

651 

79 

^s 

117 500 

108 900 

95 500 

94100 

93 700 


* Superscripts (1) and (2) refer to one-sigma (la) perturbations obtained by using the 


sensitivity-partial-derivative method and the parameter method, respectively. 
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TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION - Continued 


(b) Kp =0; Wl = w 3 = 1, = 10"^ (Set (2)) 


Sensitivity partials 

Values for sensitivity weighting factor, w s , 

of - 

and related terms 
(*) 

0 

0.1 

0.3 

1.0 

10.0 

8u/8T 

1.106 

1.084 

0.723 

0.448 

0.104 

3u/3m 

-6 322 

-5 869 

-2 267 

-103 

174 

3u/3m Q 

-84.98 

-81.22 

-43.76 

-18.88 

-3.40 


110.1 

104.4 

52.6 

22.1 

4.5 

*»g 

110.2 

104.4 

52.6 

22.1 

4.5 

dy/dT 

52.75 

52.75 

52.75 

52.75 

52.75 

8y/8m 

-96 100 

-102 100 

-108 500 

-112 100 

-126 400 

8y/8m 0 

-2 755 

-2 793 

-2 833 

-2 856 

-2 946 

< 

3 209 

3 260 

3315 

3 347 

3 474 

*g 

3 210 

3 260 

3 316 

3 347 

3 475 

8v/8T 

0.805 

0.805 

0.805 

0.805 

0.805 

8v/8m 

-1419 

-1316 

-1206 

-1 146 

-899 

8v/8m 0 

-41.74 

-41.09 

-40.40 

-40.02 

-38.46 


48.57 

47.71 

46.82 

46.35 

44.45 

to® 

48.57 

47.72 

46.82 

46.33 

44.40 

^P 

4 423 

4 336 

2 891 

1791 

417 

4>s 

1 

15 510 

14 240 

6 060 

3 760 



3 200 


* Superscripts (1) and (2) refer to one-sigma (la) perturbations obtained by using 
the sensitivity-partial-derivative method and the parameter method, respectively. 
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w s ~ 1.0 in figure 3. Part of the reason is the change in the magnitude of <p s due to 
changing only the relative values of the relative weighting factors without regard to their 
magnitude. This is clearly illustrated by the near-order-of-magnitude difference in <p s 
for set (1) and set (2) at w s = 0. (Compare tables V(a) and V(b).) The trajectories for 
the two sets are the same; <p s changes with the change from set (1) to set (2). Reducing 
the values of (p s in table V(a) by dividing by the constant ^5 510 ~ so that both 

cases have the same 0 S at w s = 0, and incorporating this constant into the weighting 
factor by multiplying each w s of table V(a) by 7.57, allows a more reasonable compari- 
son. This comparison is shown in figure 7 where the ordinate is called w s (adjusted). 
This figure shows the differences which come about due to different relative-weighting- 
factor sets. 

The control time histories for relative-weighting-factor set (2) are shown in fig- 
ure 8. Again there is a tendency toward a bang -bang type of control, as w g increases, 
which may be noted by observing that the angle difference between the nearly constant 
attitude portion of the control time history at the beginning and near the end of flight is 
about 180° for both w s = 1.0 and w s = 10.0. Figure 9 shows that in comparison with 
the set (2) trajectory for w s = 0, the other set (2) trajectories are generally less steep. 
The opposite result is shown for set (1) in figure 5; in comparison with the set (1) trajec- 
tory for w g = 0, the other set (1) trajectories are more steep. These results indicate 
the importance of the relative weighting factors in shaping the trajectories. 

Data for a nonzero value of Kp>, Kp> = 0.001, and relative -weighting -factor set (2) 
are shown in table V(c) and figures 10 and 11. Data for w g = 0.1 are shown in table V(c) 
but not in figures 10 and 11 because the control time histories and trajectories essen- 
tially coincide with those for w s = 0. In table V(c) it may be observed that while consis- 
tent data trends are shown in the first three w s columns, the data in the column for 
w g = 10.0 are not consistent. The consistent data trends in the first three columns are 
very similar to those in table V(b). The inconsistency of the last column may be explained 
by examination of figures 10 and 11 where the radically different character of the control 
time history and trajectory for w s = 10.0 may be seen. This result for w s = 10.0 was 
so unusual that it was believed necessary to verify the answer. Accordingly, verification 
was obtained by iterating to the same result (essentially) from an additional two different 
nominal trajectories. These results therefore appear correct. It is believed that the 
data inconsistency results from the different character of the trajectory - that is, the 
bending back of the trajectory as seen in figure 11. With the exception of this w g = 10.0 
case, the trajectories are less steep with increasing wg just as they were in figure 9 
for a similar case with zero drag. 

For all of the cases discussed herein, good agreement was observed (see 
tables V(a), V(b), and V(c)) between the one-sigma perturbations computed by using 
sensitivity partials and those computed by using the parameter -perturbation method. 
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TABLE V.- SENSITIVITY PARTIALS AND RELATED INFORMATION - Concluded 


(c) K d = 0.001; Wl = w 3 = 1, w 2 = 10 -4 (Set (2)) 


Sensitivity partials 
and related terms 
(*) 

Values for sensitivity weighting factor, w 

S’ of - 

0 

0.1 

1.0 

10.0 

9u/9T 

25.35 

25.24 

24.35 

12.07 

9u/3m 

-756.9 

-747.1 

-668.1 

-002.8 

8u/8m 0 

-8.11 

-7.99 

-7.04 

-1.27 

8u/8K d 

-687 000 

-688 900 

-697 500 

-359 900 


17.37 

17.28 

16.52 

6.61 

< 

17.37 

17.28 

16.52 

6.61 

8y/8T 

31.94 

31.96 

32.11 

31.81 

8y/9m 

-44 640 

-44 840 

-46 370 

-44 800 

8y/8m 0 

-1438 

-1440 

-1456 

-1428 

9y/9K D 

-14190 000 

-14 200 000 

-14 220 000 

-14 350 000 

<1 

1730 

1732 

1751 

1712 


1729 

1732 

1751 

1712 

8v/8T 

16.18 

16.13 

15.81 

17.01 

9v/9m 

-7.8 

4.6 

101.8 

34.5 

8v/9m 0 

-3.80 

-3.68 

-2.81 

-3.52 

8v/8Kj) 

-278 700 

-280 400 

-293 500 

-313 700 

< 

7.358 

7.309 

7.102 

7.766 


7.357 

7.311 

7.100 

7.766 


1824 

1824 

1811 

911 

4> s 

655 

652 

630 

403 


Superscripts (1) and (2) refer to one -sigma (Ict) perturbations obtained by using 
the sensitivity -partial -derivative method and the parameter method, respectively. 


23 



CONCLUDING REMARKS 


It has been shown how a sensitivity measure, composed of the weighted sum of the 
variances of the performance index and terminal constraints, may be added to the perfor- 
mance index of a stochastic optimal control problem to achieve a reduction in performance 
and constraint sensitivity due to parameter perturbations. 

It was necessary to assume that the parameters of the system remained fixed during 
system operation and that the stochastic nature of the problem came about because of 
inexact knowledge of these fixed values. 

It was shown how this technique increased the dimensionability of the optimization 
problem and introduced weighting factors or constants for use as design parameters. 

These weighting factors permit different relative importance to be attached to different 
types of sensitivity, such as position relative to velocity, and allow for adjustment of 
performance degradation and of sensitivity reduction. 

The feasibility of solving nonlinear multiparameter problems by using this technique 
was illustrated by solving an example rocket -trajectory problem through application of a 
steepest-descent algorithm. The example solutions also served to illustrate the tradeoffs 
made possible by changes in the weighting factors. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., April 14, 1971. 
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APPENDIX A 


STEEPEST -DESCENT DERIVATION 

The problem is to find the control time history u(t), t Q = t ^ tf, which minimizes 
4> = 4>(x (tf)) (Al) 

subject to the differential constraints 

x = f(x,u,t) ; x(t 0 ) =x Q (given) (A2) 

and the terminal constraints 


V((x(t f )) = L 


(A3) 


The solution is obtained iteratively. Choose a reasonable time history u*(t) and 
obtain x*(t), a solution to equation (A2). This solution, in general, is such that neither 
<fi is minimum nor \)s = L. 

Linearizing about this solution gives 


x 






) 


or 

Of Of 

8x = -^(t)6x + -^(t)6u (A4) 

Multiplying equation (A4) by A(t), a variable matrix integrating factor, and integrating 
by parts yields 


A6x 


** = f tf + a|£\ 6x dt + f tf A^ 6u dt 

to J t n V dt ~ 8 h - 


(A 5) 


Now, for convenience, let A(t) be the solution of 


^ + A— = 0 
dt 9x 


(A6) 
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APPENDIX A - Continued 


with 

f+> 

diL/ 

A M-sC*) 

where 



Then equation (A5) becomes 

f - A (*o)«5(to)= A 4 6sdt 

or, letting 6^ s — 6x (tf) , 

6^ = T = A ( t o) ^ A 5 ^ dt 


(A7) 


(A8) 


(A9) 


Partition A(t) as follows: 

[x 0 T ( t) 

A(t) " [A^(t')"_ 

Then, since x(t D ) is given and 6x(t 0 ) = 0, 


6u dt 


and 


. , r fc f . t 9 i 

60 “ \ X( t> 9u 
l O - 

rtf 8f 

5 ! = \ A 4's: t i it 

— 


(A10) 


(All) 


Now for some measure of allowable control perturbation 

(dP) 2 = ^ tf 6u T W 6u dt ; W - W T > 0 (A12) 

minimize 5 <j> and choose 6jp such that i^(x(t f )) = L will be satisfied (or more nearly 


so). Form an augmented function to be minimized 


J A = 50 + JLL 


rtf rp ~ I rp ( rtf 8f 

(dP) 2 - J 8u X W 6u dt + V 1 M - A^— 5u dt 


(A 13) 
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APPENDIX A - Concluded 


where ju and u are constant Lagrange multipliers. For to be an extremum, its 

variation with respect to 6u must be zero — that is, 


«A | - 2 M « =■ 0 


This requirement implies that 


rp 9f rp T 9f 

— - 2p Su W - v 1 Aw / — = 0 
v 3u — — ^ oh 


for all t. Solving for 6u, 


6 !1 = ^ W ' 1 (1) T ( X 0- A ^) 

In order to solve for ju and v this expression for 6u is substituted into equa- 
tions (A10) and (All) to obtain 


(A 14) 


^ dip + I xf/xjy ^1 \p(p 


2n = -{( 1^ - V^%^‘\ 0 )/|(dP) 2 - dip 7 !^- 1 d*j} 


l/2 


(A 15) 


where the sign on the radical is minus because <f> is being minimized and where 

*\ 


ftf T3f _i/9f\ T 

I , , B \ XA ) 1 — w 1 — X . dt 
90 A 9 8u \9u / 9 

ptf 9f i/9f\ T 

I, , . b \ Aw, — W -1 (— Xa. dt 
A ^ 8u \8u ) 9 

ftf 9f _ 1 / af \ T x 

wl b V® 

l O — \ — / 


(A 16) 


Now let 


u(t) = u (t) + 6u(t) 


be the new control time history used to obtain a solution to equation (A2) and repeat the 
same procedure until ip - L and the gradient of the payoff -constraint surface 

(AID 

are sufficiently close to zero. 
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Figure 2.- A comparison of the initial guessed thrust -attitude 
control time history with the steepest-descent generated 
optimal thrust -attitude control time history. Kp> = 0; w g = 0. 
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10 


its and sensitivity and performance indices 
tor. Kq = 0; wj = W3 = 1, W2 = 10 "2. 


Thrust-attitude control angle, 0, deg 



Time from lift-off, t, sec 

Figure 4.- Comparison of optimal control time histories for different sensitivity 
weighting factors Kr> = 0; wi = wo = 1, W 2 = 10”^. 






Percent of value at ws = 



o o 






Thrust-attitude control angle, 6 ; deg 



Figure 10.- Comparison of optimal control time histories for different 
sensitivity weighting factors. Kq = 0.001; wj = W3 = 1, W2 = 10 - ^. 


39 



National Aeronautics and Space Administration 
Washington, D. C. 20546 


OFFICIAL business 

PENALTY FOR PRIVATE USE $300 


FIRST CLASS MAIL 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 


11U 001 44 51 30S 71166 00903 

AIR FORCE WEAPONS LABORATORY /WL0L/ 
KIRTLAND AFB» NEW MEXICO 87117 


ATT E. L00 BOWMAN* CHIEF»TECH. LIBRARY 


POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereQf.” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. . 

CONTRACTOR REPORTS: Scientific and 
technical information generated u-nder a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs., 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 


